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Abstract. 

Recently, we used the Sine collocation method with the double exponential transformation to compute 
eigenvalues for singular Sturm-Liouville problems. In this work, we show that the computation complexity 
of the eigenvalues of such a differential eigenvalue problem can be considerably reduced when its operator 
commutes with the parity operator. In this case, the matrices resulting from the Sine collocation method are 
centrosymmetric. Utilizing well known properties of centrosymmetric matrices, we transform the problem 
of solving one large eigensystem into solving two smaller eigensystems. We show that only of all 

components need to be computed and stored in order to obtain all eigenvalues, where 2N + 1 corresponds to 
the dimension of the eigensystem. We applied our result to the Schrodinger equation with the anharmonic 
potential and the numerical results section clearly illustrates the substantial gain in efficiency and accuracy 
when using the proposed algorithm. 
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1 Introduction 


In science and engineering, differential eigenvalue problems occur abundantly. Differential eigenvalue prob¬ 
lems can arise when partial differential equations are solved using the method of separation of variables. 
Consequently, they also play an important role in Sturm-Liouville (SL) differential eigenvalue problems [Tj. 
For example, the solution of the wave equation can be expressed as the sum of standing waves. The fre¬ 
quencies of these standing waves are precisely the eigenvalues of its corresponding Sturm-Liouville problem. 
Similarly, in quantum mechanics, the energy eigenvalues associated with a Hamiltonian operator are mod¬ 
elled using the time-independent Schrodinger equation which is in fact a special case of a Sturm-Liouville 
differential eigenvalue problem. 

Recently, collocation and spectral methods have shown great promise for solving singular Sturm-Liouville 
differential eigenvalue problems [^|3]. More specifically, the Sine collocation method (SCM) [3H5] has been 
shown to yield exponential convergence. During the last three decades the SCM has been used extensively 
to solve many problems in numerical analysis. The applications include numerical integration, linear and 
non-linear ordinary differential equations, partial differential equations, interpolation and approximations to 
functions [ail]. The SCM applied to Sturm-Liouville problems consists of expanding the solution of a SL 
problem using a basis of Sine functions. By evaluating the resulting approximation at the Sine collocation 
points separated by a fixed mesh size h, one obtains a matrix eigenvalue problem or generalized matrix 
eigenvalue problem for which the eigenvalues are approximations to the eigenvalues of the SL operator. 
In 13, we used the double exponential Sine collocation method (DESCM) to compute the eigenvalues of 
singular Sturm-Liouville boundary value problems. The DESCM leads to a generalized eigenvalue problem 
where the matrices are symmetric and positive-definite. In addition, we demonstrate that the convergence of 
the DESCM is of the rate O ^ io^(jv)^ for some k > 0 as ^ oo, where 2iV-l- 1 is the dimension 

of the resulting generalized eigenvalue system. The DESCM was also applied successfully to the Schrodinger 
equation with the anharmoiiic oscillators nnj. 

In the present contribution, we show how a parity symmetry of the Sturm-Liouville operator can be conserved 
and exploited when converting our differential eigenvalue problem into a matrix eigenvalue problem. Indeed, 
given certain parity assumptions, the matrices resulting from the DESINC method are not only symmetric 
and positive definite; they are also centrosymmetric. The study of centrosymmetry has a long history 
dlHIH]. However, the last two decades has stemmed much research focused on the properties and applications 
of centrosymmetric matrices ranging from iterative methods for solving linear equations to least-squares 
problems to inverse eigenvalue problems [191133] . 

Using the eigenspectrum properties of symmetric centrosymmetric matrices presented in we apply 
the DESCM algorithm to Sturm-Liouville eigenvalue problems and demonstrate that solving the resulting 
generalized eigensystem of dimension (2iV -I- 1) x {2N -I- 1) is equivalent to solving the two smaller eigensystems 
of dimension N x N and (A^-|- 1) x (A^-|- 1). Moreover, we also demonstrate that only of all components 
need to be stored at every iteration in order to obtain all generalized eigenvalues. To illustrate the gain 
in efficiency obtained by this method, we apply the DESCM method to the time independent Schrodinger 
equation with an anharmonic potential. Furthermore, it is worth mentioning that research concerning inverse 
eigenvalue problems where the matrices are assumed centrosymmetric has been the subject of much research 
recently |^l28j . Consequently, the combination of these results and our findings could lead to a general 
approach for solving inverse Sturm-Liouville problems. 

All calculations are performed using the programming language Julia and all the codes are available upon 
request. 
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2 Definitions and basic properties 


The sine function valid for all z G C is defined by the following expression: 

sin(7rz) 


sinc(z) = 


TTZ 

1 


for z^O 
for z = 0. 


( 1 ) 


For j G Z and h a positive number, we define the Sine function S(j, h)(x) by: 

S{j,h){x) = sine a: G C. (2) 

The Sine function defined in @ form an interpolatory set of functions with the discrete orthogonality 
property: 

S{j,h){kh) = 6j^k for (3) 

where Sj^k is the Kronecker delta function. 

Definition 2.1. Given any function v defined everywhere on the real line and any h > 0, the symmetric 
truncated Sine expansion of v is defined by the following series: 

N 

CN{v,h){x)= ^ Vj^hS{j,h){x), (4) 

j=-N 


where Vj^h = v{jh). 

The Sturm-Liouville (SL) equation in Liouville form is defined as follows: 

Lu{x) = —u"{x) + q{x)u{x) = Xp{x)u{x) 

a < x < b u{a) = u{h) = 0, (5) 

where —oo < a < b < oo. Moreover, we assume that the function q{x) is non-negative and the weight 
function p(x) is positive. The values A are known as the eigenvalues of the SL equation. 

In |M], we apply the DESCM to obtain an approximation to the eigenvalues A of equation ([S]). We initially 
applied Eggert et al.’s transformation to equation ([5]) since it was shown that the proposed change of variable 
results in a symmetric discretized system when using the Sine collocation method [35| . The proposed change 
of variable is of the form [351 Defintion 2.1]: 

v{x) = o(j){x) u{x) = (6) 

where (j)~^(x) is a conformal map of a simply connected domain in the complex plane with boundary points 
a b such that = —oo and (j)~^{b) = oo. 

Applying the change of variable ([S]) into equation (O, one obtains |35| : 

Cv{x) = —v"(x) -I- q(x)v(x) = Xp((f(x))((p'(x))^v(x) with lim v(x) = 0, (7) 

|ai|—>-oo 

where: 

q(x) = -x/<f'(x) ^ (8) 

To implement the double exponential transformation, we use a conformal mapping (f(x) such that the solution 
to equation © decays double exponentially. In other words, we need to find a function 4>(x) such that: 

|z;(a:)| < Aexp(-Bexp( 7 |a:|)), (9) 
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for some positive constants A, B, 7 . Examples of such mappings are given in [M1I36] . 

Applying the SCM method, we obtain the following generalized eigenvalue problem: 

£Cm(i’,^) = Av = /iD^v =4> (A —/iD^)v = 0, (10) 


where the vectors v and CmIv, h) are given by: 

v= {v-N,h, - ■ ■ ,VN,hf and CN{v,h) = {CN{v,h){-Nh), ... ,CN{v,h){Nh))'^ , 


( 11 ) 


and p are approximations of the eigenvalues A of equation ([7]). For more details on the application of the 
SCM, we refer the readers to [H]- 

As in we let denote the Sine differentiation matrix with unit mesh size: 


liI '■><“=> 


( 12 ) 


x—kh 


The entries Aj^k of the {2N + 1) x {2N + 1) matrix A are then given by: 

^hk = -^5fj:+q{kh)5fl with -N<j,k<N, (13) 

and the entries Djf. of the {2N + 1) x {2N + 1) diagonal matrix are given by: 

Dl„ = {^'{khjfp{^{khj)6l°l with -N<j,k<N. (14) 


As previously mentioned, Eggert et al.’s transformation leads to the matrices A and to be symmetric 
and positive definite. However, as will be illustrated in the next section, given certain parity assumptions, 
these matrices yield even more symmetry. 


3 Centrosymmetric properties of the matrices A and 


In this section, we present some properties of the matrix A and that will be beneficial in the computation 
of their eigenvalues. The matrices A and are symmetric positive definite matrices when equation © is 
discretized using the Sine collocation method. Additionally, given certain parity assumptions on the functions 
q{x), 4){x) and p{x) in equation ([T]), the matrices A and will also be centrosymmetric. 

Definition 3.1. ]3iA Section 5.10] Let J denote the parity operator defined by: 

Jfix) = f{-x), (15) 

where fix) is a well defined function being acted upon by J. 

Definition 3.2. An operator B is said to commute with parity operator if it satisfies the following relation: 

BJfix) = JBfix). (16) 

Equivalently, we can say that the the commutator between B and J is zero, that is: 


[B, J]=BJ- JB = 0. 


(17) 


Definition 3.3. fSOl Definition 5] An exchange matrix denoted by 3 is a square matrix with ones along the 
anti-diagonal and zeros everywhere else: 


J = 


/o 

Vi 



(18) 
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Definition 3.4. fc)9l Definition 2] Let B &e a matrix of dimension (2N + 1) x {2N + 1) with components 
Bj^k for —N < j,k < N. B is centra symmetric if and only i/ B satisfies the following property: 

BJ = JB, (19) 

where J is an exchange matrix of dimension (2N + 1) x (2N + 1). Writing equation (1191) in a component 
form, we have the following relation: 

B-j-k = Bj^k for - N < j, k < N. (20) 


We now present the following Theorem establishing the connection between symmetries of the Sturm- 
Liouville operator and its resulting matrix approximation. 

Theorem 3.5. Let C denote the operator of the transformed Sturm-Liouville problem in equation 

picfixm'ixw 

If the commutator = 0, where J is the parity operator, then the matrices A and defined by 

equations (USD and m resulting from the DESCM are centra symmetric. 


Proof The commutator [C,U] = 0 if and only if q(x) and p{x) are even functions and 4>{x) is an odd 
function. 


If <j){x) is an odd function, then (j)' {x) is even, (j)''{x) is odd and 4)"'{x) is even. From this and equation ([5]), 
it follows that q{x) is even. 

In order to show that the resulting matrices A and are centrosymmetric, we demonstrate that both these 
matrices satisfy equation (HOD- Before doing so, it is important to notice that the /*** Sine differentiation 
matrices defined in equation (1121) have the following symmetric properties: 



S{-j, h){x) 


x— — kh 


Sj‘l if I is even 
if I is odd. 


( 22 ) 


Hence, the /*** Sine differentiation matrices are centrosymmetric if I is even. It is worth noting that when I is 
odd, the Sine differentiation matrices are skew-centrosymmetric m- Consequently, investigating the form 
for the components of the matrix A in equation (1131) . we obtain: 

+ q{-kh) 

= (23) 

Similarly, investigating the form for the components of the matrix in equation (I14L we obtain: 

= i^'{-kh)fpm-kh))S^°l_, 

= {<!}'ikh)f p{(j){kh)) Sfl 

= Dl,. (24) 

Both matrices A and satisfy equation (EPl) . From this it follows that A and are centrosymmetric. 

Theorem [3^ illustrates that Sine basis functions preserve the parity property of the Sturm-Liouville operator 
when discretized. Hence, when the matrices A and are symmetric centrosymmetric positive definite 
matrices, we can utilize these symmetries when solving for their generalized eigenvalues. In m. Cantoni 
et al. proved several properties of symmetric centrosymmetric matrices. In the following, we will utilize 
some of these properties to facilitate our task of obtaining approximations to the generalized eigenvalues of 
the matrices A and D^. The following lemma will demonstrate the internal block structure of symmetric 
centrosymmetric matrices. 
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Lemma 3.6. Lemma 2] If H is a square symmetric centra symmetric matrix of dimension {2N + 1) x 
(2iV + 1), then H can be written as: 



S 

X 

CT ■ 

H = 

T 

X 

h 

xTj 


C 

Jx 

JSJ 


(25) 


where S, C are matrices of size N x N, J is the exchange matrix of size N x N, x. is a column vector of 
length N and h is a scalar. In addition, S'^ = S and C'^ = JCJ. 


The next lemma simplifies the calculation needed to solve for these eigenvalues. 

Lemma 3.7. ]12[ Lemma 3] Let H be a square symmetric centrosymmetric matrix as defined in lemma [Kd\ 
and let 'V be a square matrix of dimension (27V + 1) x {2N + 1) defined by: 


V = 


S-JC 0 0 

0 h 1/2 x'^ 

0 V2x S + JC 


(26) 


then H and V are orthogonally similar. That is, the matrices H and V have the same Jordan normal form 
and thus the same eigenvalue spectrum. 


Canton! et al. use Lemmas 13.61 and 1X71 to prove the following Theorem concerning a standard eigenvalue 
problem where the matrix is centrosymmetric. 


Theorem 3.8. J121 Theorem 2] Let H he a square symmetric centrosymmetric matrix as defined in Lemma 
[M then solving the eigenvalue problem det(H — AI) = 0 is equivalent to solving the two smaller eigenvalue 
problems: 


det(S — JC — AI) = 0 and det 


h \f2 x"^ 
y2x S + JC 


-A 


1 0 

0 I 


= 0 . 


(27) 


Since our problem consists of solving a generalized eigenvalue problem where one matrix is a full symmetric 
centrosymmetric and the other is a diagonal centrosymmetric matrix, we propose the following Theorem. 

Theorem 3.9. Let H and W be square symmetric centrosymmetric matrices of the same size, such that: 

(28) 


H = 

S 

T 

X 

X 

h 

1 

and 

W = 

diag (w) 

0 

0 

w 

0 

0 


C 

Jx 

JSJ 



0 

0 

Jdiag (w)J 


then solving the generalized eigenvalue problem det(H — AW) = 0 is equivalent to solving the two smaller 
generalized eigenvalue problems: 


det(S — JC — Adiag (w)) = 0 and det 


h \f2 x'^ 
72 X S + JC 


-A 


w 0 

0 diag (w) 


= 0 . 


(29) 


Proof This proof relies on the unitary transformation matrix presented in Lemma 3]: 

1 


K = 


72 


I 0 -J 

0 72 0 

I 0 J 


(30) 


where I is the identity matrix and J is the exchange matrix. 
It is easy to verify that: 

KHK'^ = V, 

where V is the matrix in lemma 1X71 


(31) 
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This result is analogous for the matrix W with a change in notation. Hence: 


0=det(H-AW) 

= det(K) det(H - AW) det(K'^) 
= det(KHK'^ - AKWK'^) 
S-JC 0 

= det I 0 h 

0 X 


0 

72 xT 

S +JC 


= det(S — JC — Adiag (w)) det 


-A 


diag (w) 

0 

0 


h 

72 X 


72 xT 

S +JC 


w 

0 


-A 


0 
0 

diag (w) 

w 0 

0 diag (w) 


(32) 


from which the result follows. 

Theorem 13.91 is very useful when N is large since it is less costly to solve two symmetric generalized eigen- 
systems of dimensions N x N and {N -I- 1) x {N + 1) rather than one symmetric eigensystem of dimension 
{2N -I- 1) X (2N -I- 1). Additionally, Lemma [3^ also has large ramifications when it comes to saving storage 
space. As is discussed in |3Il, the Sine differentiation matrices are symmetric toeplitz matrices. Therefore, 
for a symmetric toeplitz matrix of dimension {2N -I- 1) x {2N + 1), only 2N -|- 1 elements need to be stored. 
Investigating the definition of the matrix A in equation (fT^ . we can see that A is defined as the sum of a 
symmetric toeplitz matrix and a diagonal matrix. Moreover, from Lemma 13.61 and Theorem 13.91 using only 
the antidiagonal and anti-upper triangular half of matrix C, the vector x, the scalar h, the diagonal and 
lower triangular half of the matrix S, the vector w and the scalar w, we can create all the elements needed 
to solve for the generalized eigenvalues of the matrices A and D^. Hence, the ratio of elements needed to 
be computed and stored at each iteration N in order to solve for these eigenvalues is given by: 


Proportion of Entries Needed 


{2N + N+1) + {N+1) 
{2N+iy + {2N+1) 


1 

N+1' 


(33) 


Thus, only entries need to be generated and stored at every iteration to obtain all of the 

generalized eigenvalues. 

In the following section, we will illustrate the gain in efficiency of these results by applying the DESCM to 
the Schrodinger equation with an anharmonic oscillator. 


4 The anharmonic oscillator 

The time independent Schrodinger equation given by: 

'Hip[x) = E'tp{x) with lim '!/;(x) = 0. (34) 

|ai|^oo 


In equation (1341) . the Hamiltonian is given by the following linear operator: 

d^ 

where V{x) is the potential energy function and E is the energy eigenvalue of the hamiltonian operator T-L. 
In our case, we are treating the anharmonic oscillator potential V(x) defined by: 

m 

V{x) = with Cm > 0 and m G N\{1}. (35) 

i=l 

In |10| . we successfully applied the DESCM to time independent Schrodinger equation with an anharmonic 
potential. As we can see, the time independent Schrodinger equation (1341) is a special case of a Sturm- 
Liouville equation with q{x) = V{x) and p{x) = 1. Applying Eggert et al.’s transformation and the DESCM 
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with (/)(x) = sinh(a:), we arrived at the following generalized eigenvalue problem: 


det(A-fD2) = 0, (36) 

where £ are approximations of the energy eigenvalues E. 

The matrices A and defined by equation (fT51) and (I14|) are given by: 

Aj,k = - Sfl + V{kh)5fl with -N<j,k<N, (37) 

where: 

13 ™ . 

V(x) = ^ ~ ^ sech^(a:) + cosh^(x) ^ q sinh^*(a:), (38) 

i=l 

and: 

= cosh^(fc/i) (5^*2 with — N < j,k < N. (39) 

Since the anharmonic potential V (x) defined in (1351) is an even function, the function p{x) = 1 is an even 
function and the conformal map (/)(a:) = sinh(a;) is an odd function, we know that Theorem 13.51 applies. 
Hence, the matrices A and are symmetric centrosymmetric. 


5 Numerical Discussion 


In this section, we test the computational efficiency of the results obtained in Theorem l3.9l All calculations 
are performed using the programming language Julia in double precision. The eigenvalue solvers in Julia use 
the linear algebra package LA PACK. 

In [40], Chaudhuri et al. presented several potentials with known analytic solutions for energy levels calcu¬ 
lated using supersymmetric quantum mechanics, namely: 

Vi{x) = x^ — Ax'^ + => 

V2{x) = Ax^ — dx'^ + x^ => 

Vslx) = (105/64)x^ — (43/8)x^-I-a:® — a;®-I-=> 

Vilx) = (169/64)a:^ — (59/8)a:^-t — a®-I-=> 

Figure [T] presents the absolute error between our approximation and the exact values given in (1401) . The 
absolute error is defined by: 

Absolute error = \£i{N) — Exact value] for Z = 0,1. (41) 


Eo = -2 
El = -9 
Eo = 3/8 
El = 9/8. 


(40) 


The optimal mesh size obtained in |10j : 


h = 


W 


2'^7r^(m+l)N \ 


(m + 1)N 


(42) 


where W (x) is the Lambert W function, is used in the calculation. 

As can be seen from Figure [H using the centrosymmetric property improves the convergence rate of the 
DESCM significantly. 


6 Conclusion 

Sturm-Liouville eigenvalue problems are abundant in scientific and engineering problems. In certain appli¬ 
cations, these problems possess a symmetry structure which results in the Sturm-Liouville operator to be 




commutative with the parity operator. As was proven in Theorem 13.51 applying the DESCM will preserve 
this symmetry and results in a generalized eigenvalue problem where the matrices are symmetric centrosym- 
metric. The centrosymmetric property leads to a substantial reduction in the computational cost when 
computing the eigenvalues by splitting the original eigenvalue problem of dimension (27V + 1) x (27V + 1) 
into two smaller generalized eigensystems of dimension TV x TV and (TV + 1) x (TV + 1). Moreover, due to 

the internal block structure of the matrices obtained using the DESCM, we have shown that only ——- 

of all entries need to be computed and stored at every iteration in order to find all of their eigenvalues. 
Numerical results are presented for the time independent Schrddinger equation (1341) with an anharmonic 
oscillator potential (1551) . Four exact potentials with known eigenvalues are tested and the results clearly 
demonstrated the reduction in complexity and increase in convergence. 


7 Tables and Figures 


Absolute Error for energy level 0 


Absolute Error for energy level 1 



(a) 


Absolute Error for energy level 0 




(b) 


Absolute Error for energy level 1 



(c) 


(d) 


Figure 1: Absolute error for the potentials Vi{x) for i = 1, 2, 3,4 given by (1501) with 4){x) = sinh(x). 

(a) Vi{x) = x^ — 4x‘* + x^ with exact eigenvalue Eq = —2. (b) V 2 (x) = ix^ — 6a;^ + x® with exact 

eigenvalue Ei = —9. (c) V 3 (x) = (105/64)x^ — (43/8)x^ + x® — x® + x^° with exact eigenvalue Eg = 3/8. 
(d) V 4 (x) = (169/64)x^ — (59/8)x‘^ + x® — x® + x^° with exact eigenvalue Ei = 9/8. 
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